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We perform a variational quantum Monte Carlo simulation of an interacting Fermi gas confined 
in a three dimensional harmonic potential. This gas is considered as the precursor system from 
which a molecular bosonic gas is formed. Based on the results of two-body calculations for trapped 
atoms, we propose a family of variational many-body wave functions that takes into account the 
qualitative different nature of the BCS and EEC regimes as a function of the scattering length. 
Energies, densities and correlation functions are calculated and compared with previous results for 
homogeneous gases. Universality tests at the unitarity limit are performed including the verification 
of the virial relation and the evaluation of the /3 parameter. 
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The crossover from a low interacting attractive 
Bardeen-Cooper-Schieffer (BCS) gas to a molecular Bose- 
Einstein condensate (BEC) has been realized experimen- 
tally using a mixture of ultracold Fermi atoms in two 
hyperfine states [l], 0, S Hi- In dilute Fermi gases, 
the atomic interactions have a range much smaller than 
the interparticle separation. Nevertheless, under proper 
conditions, a magnetic field can be used to tune the at- 
tractive potential so that the energy of a pair of scat- 
tering atoms is close to that of a molecular bound state 
(Feshbach resonance). In experiments where the reso- 
nance is broad, it can be represented by a single-channel 
model in which the s-wave scattering length a determines 
the general features of the ultracold atomic gas. At low 
enough temperatures, atoms in different hyperfine states 
pair into bound molecules for positive values of a, form- 
ing a molecular BEC that can be adiabatically converted 
into a degenerate Fermi gas by shifting a to negative val- 
ues. At resonance (|a| oo), the gas acquires universal 
properties, i. e., they are independent of any feature of 
the atomic potentials 0, @, 01 ■ This is the so called uni- 
tarity limit. 

The theoretical description of the BEC-BCS crossover 
is complicated because there is no ad hoc single param- 
eter to control the relevant features in both regimes. A 
reasonable alternative to infer properties of the system 
at unitarity has been to consider an homogeneous Fermi 
gas, and to assume the universality hypothesis according 
to which the only dominant length is the interparticle 
separation. Then, the thermodynamic potentials have a 
universal form specified by only few universal numbers 
0. For instance, the interaction energy is proportional 
to the Fermi energy via a universal constant /3. 

Quantum Monte Carlo (QMC) techniques support ah 
initio methods to theoretically test the universality hy- 
pothesis. They can be used to approximately solve the 
many-body Schrodinger equation for a given model of the 
interaction potential. In previous studies the BEC and 
BCS regimes have been explored using a fixed-node QMC 
technique, which is rigorous but also computationally de- 



manding. In particular, the value of /3 has been predicted 
considering up to 66 particles [1, 0]. Nevertheless, these 
QMC calculations have been performed for homogeneous 
gases despite the fact that experimentally the atoms are 
confined by an external field and therefore the inhomo- 
geneity is intrinsic to the problem. Comparison of QMC 
results with experiments are then based on a mapping of 
the trapped system onto a corresponding homogeneous 
problem with a local value of the Fermi energy ei?(x) (lo- 
cal density approximation). This reasoning also allows 
to study the universality hypothesis consequences for a 
trapped dilute Fermi gas in the crossover 0, [l^, [ll| . 



Here we study the behavior of a balanced mixture of 
2N interacting trapped Fermi atoms, using a variational 
quantum Monte Carlo (VQMC) technique, which allows 
to deal with large number of particles but strongly de- 
pends on the choice of the variational wave function. We 
compute the many-body ground state of the inhomoge- 
neous gas from the BCS to the BEC regime for a given 
short range interaction potential. For simplicity we shall 
consider an isotropic harmonic potential. The values of 
N are the highest reported in QMC calculations for this 
kind of atomic systems. Emphasis will be made on the 
energies, densities and two point correlation functions. 
At the unitarity limit, (3 will be directly evaluated and 
the validity of the virial theorem verified 



To describe the many-body system, let us first consider 
the attractive two-body interaction potential V{rij) — 

— Vqc"^''^' "'^3 of range 6/2, to model the effective in- 
teraction between atoms in different hyperfine states. 
The corresponding two-body problem for a particle of 
reduced mass m/2, [fP /m -\- V]if = E(.p, has analytical 
s-wave solutions |12| ip{rij) = u{rij)/rij both in the con- 



tinuum {u{y) = ciJ 
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in the bound region {u{y) ~ c+J. 
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Bessel function of order v. The scattering length is 
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with A'o the second kind Bessel function of zero order, 
and C the Euler constant. This scattering length di- 
verges whenever Jo{by/Vom/h) = 0. If 2:„, n = 0, 1, 2, ... 
are the zeros of this Bessel function in increasing or- 
der, the potential V{rij) admits just n-bound states 
for z„ < by/Vom/h < z„+i. The discrete eigenvalues 
are determined by the boundary condition at = 0, 



J^^ I^I^ ^^(6-\/^™/^) ~ 0. Taking into account the pres- 
ence of a spherical harmonic trap the Schrodinger equa- 
tion, [p'^/m + TO(jj^r^/4 + V{r)]ip{r) = E(p{r), can be 
numerically solved and some interesting results are ob- 
tained when the range of the potential b/2 << y^h/muj 
[l3| . For instance: (i) The ground state energy Eq for 
by^Vom/h — zq (which implies \a\ — > oo and a zero- 
energy resonance in the homogeneous problem) becomes 
Eq ^ O.btki! as 5 ^ (e.g. for b = 0.03y/h/muj, 
Eq ~ 0.510655588?ki;) while the s-wave excited states 
have energies ^ (2rie -t- 0.5)?iw; (ii) the ground state 
for zq < b\/Vom/h < z\ can be well approximated by 
u{y{r))exp{—mr^ /2huj)/r. Result (i) is consistent with 
the analytical solution of the two-body problem in the 
presence of a contact interaction with a coupling con- 
stant determined by the scattering length [l^ . 

To address the many-body system wc rely on the 
VQMC method. A trial wave function (j>T is as- 
sumed and, from initial stochastically generated events, a 
Metropolis Monte Carlo algorithm samples the distribu- 
tion |0Tp/ / |0Tp for variations of the atoms positions. 
After a thermalization process, the energy is evaluated at 
each step. The average energy converges to its expecta- 
tion value, provided enough points are taken to sample. 
The parameters used to characterize 4>t are then var- 
ied to search a minimum of the energy. The strength of 
this method relies on the proper choice of the trial wave 
function. Here, the trial wave function in the region of 
negative scattering length, < b^/Vom/h < zq, has the 
Jastrow-Slater form 



<S>fg{V^x) (2) 



where A and to' are variational parameters, and ^fg{x) is 
the Fermi gas wave function given by a product of Slater 
determinants (one for each hyperfine state) describing 
a noninteracting system of harmonically trapped atoms. 
This variational wave function has the advantage of being 
exact when no interactions between hyperfine states are 
allowed (A = 0, w' = a; ) and does not require an explicit 
introduction of a healing distance [1, ■ It is insp ired on 
previous calculations for the nuclear matter [1^], where 
an appropriate choice of the potential allows to explore 
dynamically the interplay of the nuclear-to-quark matter 
regime while being exact in both limits. Other forms of 
the Jastrow wave function can be found in Refs. [1, Q. 



In the region of non-negative scattering length and 
Zq < b^/Vom/h < 2i, the Jastrow-Slater wave function 
over estimates the energy. Therefore, we propose the fol- 
lowing trial wave function 

■i>x{x)=AMhl')-MN,N') (3) 

where A is the antisymmctrizcr operator and 

Mhj) - ii(y(r,,))e("'"l^^-l'/^''")e(-'™l^"l'""V^u (4) 

with fij and Rij the relative and center of mass coordi- 
nates associated to rj and rj ; A is the variational param- 
eter. The structure of this trial wave function is a vari- 
ational extension for the inhomogeneous gas from that 
proposed in Ref. . 

Given 2A^ particles, a fixed value of b and a scattering 
length a, we determine the energy for a set of values of the 
variational parameters to pick up the optimal. Each run 
used about 10'^ steps for thermalization and about lO** 
more to take data. In addition, we estimate the depen- 
dence on the initial conditions that might be not erased 
in the thermalization process. The quoted errors take 
into account all the above factors. 

The results obtained for the optimal energy per parti- 
cle as a function of the scattering length a are shown in 
Fig. 1. The values are normalized as described in its cap- 
tion. By construction these normalized energies take the 
asymptotic value of zero (one) for small positive (nega- 
tive) values of a. The fact that these energies fall in a 
curve almost independent of the value of N reflects that 
the Fermi wave number kp, associated to the ideal Fermi 
gas energy E'/fg = h^kp/2m « {6N)^/'^huj, defines the 
proper distance scale to measure a. It also means that 
the same curve can be expected for larger values of 
and supports a universal behavior in the crossover. 

Notice that in the BCS side of the crossover the inter- 
action energy is very small compared to Ejpc, in accor- 
dance with previous QMC calculations in homogeneous 
gases d, Q. In fact, we have checked that for —kpa < 1 
this variational energy coincides with the result obtained 
using an effective contact interaction At unitarity 

and for a > 0, the role of the interaction energy becomes 
more relevant. It is found that in the extreme BEG limit 
kpa — > 0"*" there is a strong competition between Pauli 
blocking effects and the attractive interaction between 
the atoms in different hyperfine states. This leads to nu- 
merical precision problems that could indicate that this 
highly correlated system should be described beyond the 
simple scheme of almost non-interacting molecules. 

The (3 parameter for a trapped gas is usually evaluated 
according to the following reasoning 0, . Within the 
universality hypothesis, for a trapped gas at unitarity, 
the equation of state would be (1 -I- /3)ep + Uho = MOj 
with ep{x) the local Fermi energy, Uho the trapping 
potential, and the global chemical potential. This 
equation is equivalent to that of a trapped noninteracting 
gas of particles with an effective mass m/(l -I- (3) [3 so 
that the effective chemical potential is also simply scaled 
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to Ho = EjpgV^ + The total energy at T = 0, which 
determines the energy scale is then: 



E 



IFG 



(5) 



Since we are working with few atoms it is necessary to 
take into account the exact expression for EipQ. For a 
closed shell configuration (all single particle states with 
energies below Ep = [Mf + 3/2)huj are occupied) this 
energy is given by [l7j 



Eifg/'^NHuj = (3/4)(Xf + 2), 



(6) 



instead of the large N limit, Eifq/'^N = SEp/A, which 
produces a slightly underestimated value. 

The universality hypothesis can be tested indepen- 
dently through a virial relation [ll[ resulting from me- 
chanical equilibrium conditions on the trapped unitary 
gas. According to it, for our system 



Eu/2 = 2N{muj'^R^/2). 



(7) 



In Table 1 we show the energy per particle at unitarity 
Eij/2N obtained from the VQMC calculation and the 
trap energy muj'^{R'^) for closed shell configurations with 
Ali? < 9. It can be noticed that, within error bars, the 
virial relation Eq.© is satisfied. The approximate linear 
behavior of Ejj /2N as a function oi Aip: 

Eu/2Nhuj - 0.53 ± Om{M f + 1.95± 0.06) (8) 



supports the universality relationship Eq. ([5|) with the 
universal parameter /3 = — 0.50;to;o4. 

Previous theoretical calculationspredict (3 ~ —0.326 

0, /3 0.4_Ili,[li], P 0.56 Hi], P = -0.75 

(3 = -0.492 [iirand (3 = -0.545 The first experi- 
mental measurements yielded 5 ~ —0.3 [3], (3 = — 0.49± 
0.04 [3, /3 = -0.64 ±0.15 [| and /3 = -0.68t!^:i^[3. 
More recently the value (3 = -0.54 ± 0.05[2i| and (3 = 
-0.54;° °^ H was found for ^Li and '^^ K respectively. 

The single-particle and the two-particle correlation 
functions were calculated as a function of the scattering 
length. First, wc analyze the density profile as a function 



of the distance to the center of the harmonic trap (Fig. 2). 
The numerical profile of the ideal gas reproduces satisfac- 
torily the Thomas-Fermi (TF) prediction. At unitarity, 
the density resembles more a TF than a Gaussian pro- 
file, as predicted in Ref. [1^. For the BEG regime a clear 
molecular bunching is observed around the origin, how- 
ever, due to molecular repulsion, these bosonic molecules 
are prevented from collapsing to the center of the trap. 

By analyzing the two point correlation functions, Pauli 
blocking can be observed. As expected, the radii at which 
no atoms of the same species can be found diminishes 
as the intensity of the interaction increases for a given 
range of the potential [ll]. In Fig. 3 we compare the two 
point correlation functions -fC(rJ,rj) of atoms in differ- 
ent hyperfine states in the ideal and BEG regimes. The 
enhancement of correlation for short r^- near the center 
of the trap {Rem < 0.65-\/?i/ma;) indicates molecule for- 
mation, while the increase of probability of finding pairs 
of particles se parated at relative distances of the order of 
1.3 and 2.3 ^/hj muj with center of mass radii 0.65 and 
1.09 y^h/ niLO is a manifestation of molecular condensa- 
tion. A similar analysis performed for the BGS regime 
shows a slight increase of long distance correlations of 
atoms in different states. 

Summarizing, we have performed for the first time 
direct tests of the universality hypothesis in the uni- 
tarity limit for an inhomogeneous interacting Fermi 
gas using VQMC techniques. These tests include: the 
iV— independent energy curve features as described in 
Fig. 1, the verification of virial relations for each N and 
the variational evaluation of (3 using a linear fit of the en- 
ergy per particle Ejj /2N as a function of the Fermi num- 
ber M.F- We have also shown that the optimized wave 
functions yield reasonable values compared with experi- 
mental observations not just for the energy per particle, 
but also for the mean radii at the different regimes, the 
densities and the two point correlation functions. 
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TABLE L Optimized variational parameter A, variational val- 
ues of the energy per particle Eu/2N and the trap energy 
muj{R^) for closed shells as a function of A^f at unitarity. 
The error bars in the last case take into account statistical 
sources and uncertainties in the value of A. The interaction 
potential range h/2 is varied so that kFb/2 = 0.01. 
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FIG. 1: Variational energy per particle E/2N as a func- 
tion of l/kpa. The normalization factor ejv depends on 
the number of particles and is defined as follows. The 
ideal gas energy per particle, Eifg/2N , is the asymptotic 
value of the energy per particle for small negative scatter- 
ing lengths E{a-)/2N; the asymptotic value E{a+) /2N cor- 
responds to the energy per particle for small positive scat- 
tering lengths (here a+ = 0.012 y'^/ma;) minus half the 
bound molecular energy eb{a+)/2 for this value of a (here 
ei,(a+) = 2540.22/iti^). The A'^ dependent normalization fac- 
tor ejv is the difference between E{a-) and E{a+). Finally 
€b[a)/2 = lb{a)/2 -f E{a+) /2N . In these calculations the 
range of the interactions is 6/2 = 0.0l5y' h/muj. 
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FIG. 2: Density profiles for 2N — 330 atoms. Tlie range of 
tlie potential is 6/2 — 0.015y^ h/mtu and the scattering length 
for the BEC density corresponds to l/kpa = 3.25. 
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FIG. 3: Probability difference AK that two particles with 
antiparallel spin are found separated a distance rij in the 
BEC and ideal regimes. Each curve in this figure correspond 
to different spherical radii Rem measured from the center of 
the trap. Calculations are performed at b/2 = 0.015-\/ h/mui 
and l/kpa = 3.25. 



